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I.  INTRODUCTION 


Although  much  progress  has  been  made  during  the  last  decade  in  the 
development  of  computational  techniques  for  nonlinear  analysis,  there  is 
still  a  lack  of  effective  solution  methods  for  contact  problems.  This 
is  largely  due  to  the  fact  that  the  analysis  of  contact  problems  can  be 
computationally  very  difficult,  even  for  the  simplest  geometric 
conditions  and  constitutive  relations  used.  Much  of  the  difficulty 
lies  in  that  the  boundary  conditions  of  the  bodies  under  consideration 
are  not  known  prior  to  the  analysis  but  depend  on  the  solution 
variables. 

In  two  earlier  communications  we  presented  a  solution  method  for  the 
analysis  of  general  two-dimensional  (plane  stress,  plane  strain  and 
axisymmetric)  contact  conditions  [1  2J.  The  objective  of  this  report  is  to 
present  an  algorithm  for  three-dimensional  contact  problems.  The  solution 
procedure  is  an  extension  of  the  Lagrange  multiplier/segment  algorithm, 
discussed  in  our  earlier  work,  to  three-dimensional  analysis  and  in  this 
report  we  assume  that  the  reader  is  familiar  with  the  developments  of 
References  [12]  and  the  notation  of  Reference  [3]. 

The  Lagrange  multiplier/segment  algorithm  can  be  employed  for  three- 
dimensional  large  deformation  contact  problems  with  variable  contact  areas, 
and  in  static  and  dynamic  conditions.  In  the  following  sections  we  first 
present  the  basic  equations  used  in  the  solution  algorithm  and  discuss  the 
important  numerical  details.  We  then  give  the  solution  of  various  example 
problems  to  demonstrate  the  applicability  and  limitations  of  the  algorithm 
used.  The  report  concludes  with  some  thoughts  on  future  developments  that 
should  be  pursued  to  arrive  at  improved  solution  methods  for  contact 
problems. 


II.  SOLUTION  METHOD  FOR  CONTACT  -  STATICS 

The  problem  considered  herein  deals  with  the  stress  analysis  of  two 
bodies,  contactor  and  target,  when  their  boundaries  come  into  contact  with 
each  other  under  the  action  of  external  loads  (see  Fig.  1).  The  occurrence 
of  contact  can  be  at  any  arbitrary  location  on  the  body  boundary  and  the 
basic  geometric  condition  of  contact  is  that  no  material  overlap  can  occur 
between  the  bodies.  Depending  on  the  external  loads,  large  changes  in  the 
region  of  contact  are  possible  including  relative  sliding  or  possible 
separation  after  contact.  The  developed  forces  of  contact  on  the  two  bodies 
must  be  statically  equivalent  to  each  other  and  for  each  body  the  support 
reactions  must  be  in  equilibrium  with  the  externally  applied  forces  and  the 
contact  forces. 


7 


PRESCRIBED 

FORCES 


a)  Condition  prior  to  contact 


Figure  1  Schematic  representation  of  problem  considered 
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b)  Condition  at  contact 


Figure  1  Continued 


Although  the  discussion  focuses  on  two  bodies  coining  into  contact  with 
each  other,  the  solution  method  is  also  applicable  to  the  analysis  of 
more  than  two  contacting  bodies.  In  that  case,  one  contactor  body  comes 
into  contact  with  more  than  one  target  body  and  vice  versa.  The 
calculations  for  contact  are  performed  for  each  contact  region  separately 
and  the  combined  effect  onto  the  incremental  equations  of  equilibrium  is 
obtained  by  summing  the  individual  contributions  by  the  direct  stiffness 
method. 


2.1  Discretization  of  the  Contactor  and  Target  Body  Surfaces  by 

Surface  Segments 

Both  the  contactor  and  target  surfaces  are  discretized  using  four-node 
quadrilateral  segments  (see  Fig,  lc).  Considering  the  finite  elements  which 
are  used  to  discretize  the  continuum  of  a  body,  a  generic  surface  segment 
corresponds  to  a  finite  element  face  that  lies  on  the  body  boundary. 

In  the  incremental  finite  element  solution,  the  contactor  surface  nodes 
are  considered  to  come  into  contact  with  the  target  segments.  A  contactor 
segment  is  defined  to  be  in  contact  if  all  four  nodes  belonging  to  the 
segment  are  in  contact  with  the  target  surface  (see  Figs.  lc,d). 

In  general,  the  surface  segments  are  non-flat  and  therefore  the 
following  three  geometric  assumptions  are  used  in  the  solution  procedure: 


•  For  a  generic  contactor  segment,  j,  the  normal  vector  for  the  entire 

segment  1$  approximated  by  the  exact  normal  to  the  segment  surface 

c  c 

at  (r*0,s»0),  denoted  as  n^  (see  Fig.  2a).  The  vector  n^  is 

used  in  the  calculation  of  the  total  normal  contact  force  developed 
over  segment  j  (see  Section  2.5). 


•  For  a  generic  target  segment,  j,  the  normal  vector  for  the  entire 
segment  is  approximated  by  the  exact  normal  to  the  segment  surface 

T  T 

at  (r*0,$*Q),  denoted  as  n^  (see  Fig.  2b).  The  vector  n^  is  used 

in  the  calculation  of  constraints  on  incremental  surface 
displacements  due  to  contact  (see  Section  2.2). 
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d)  Statically  equivalent  nodal  contact  forces  on  the  contactor 
and  target  surfaces 


Figure  1  Continued 


c  c 

a)  Normal  vector  to  the  contactor  segment  j  (  r£j  acts  into 

the  continuum  of  the  contactor  body) 


B 


T  T 

b)  Normal  vector  to  the  target  segment  j  (  acts  into 

the  continuum  of  the  target  body) 

Figure  2  Geometry  approx iraations  for  contactor  and  target  segments 


•  For  a  generic  target  segment,  j,  the  geometry  of  the  segment  surface 
is  approximated  by  four  triangles  which  have  one  common  vertex,  0, 
(r=0,s=0)  (see  Figs.  2b, c).  Then,  considering  a  generic  triangle  ABO, 
the  coordinates  of  any  point  P  within  the  triangle  ere  given  by 


t+At 


-P 


t+At 


iA 


+  6 


t+At 


XD  +  Y 


t+At 


*0 


(1) 


where. 


a,8,y  =  triangular  area  coordinates  of  point  P  at  time  t+At 

t+Atx.,t+AtxB  =  global  coordinates  of  target  nodes  A  and  B, 

”A  respectively  at  time  t+At 

t+Atxn  *  global  coordinates  of  vertex  0  (r-0,s=0)  at 
'u  time  t+At. 


Also,  using  the  bilinear  interpolation  functions, 


mt*0  ■  1  [t+“iA  *  tt4txB  *  t+ltJc  *  tWtxD?  <*> 


where, 


t+Atxr,  t+Atxn  a  global  coordinates  of  target  nodes  C  and  D, 
"u  respectively  at  time  t+At  . 


Substituting  Eq . (2)  into  Eq.(l),  we  obtain 

-  I  [(4a  ♦  v)t+U5A  *  <4.  *  Y)t+AtiB  *  Ytt4txc  *  Yt+4tSD]_(3) 


Equation  (3)  Is  used  in  the  calculation  of  constraints  on  incremental 
surface  displacements  due  to  contact  (see  Section  2.2). 


AREA a  Toy 


APEA=aay 


■hREA*(3gj 


TOTAL  AREA  OF 
TRIANGLE  ABO>01 


c)  Triangle  ABO  of  the  target  segment  j 


Figure  2  Continued 
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2.2  Constraints  on  Incremental  Surface  Displacements  Due  to  Contact 


After  the  two  bodies  have  come  into  contact,  the  incremental 
contactor  surface  displacements  must  be  compatible  with  the  incremental 
target  surface  displacements  so  that  the  current  conditions  of  sticking 
contact  and  sliding  contact  between  the  contactor  and  target  surfaces 
are  satisfied.  This  compatibility  of  surface  displacements  is  only 
enforced  at  the  discrete  locations  corresponding  to  the  contactor  nodes. 

As  a  result,  in  an  equilibrium  configuration,  the  contactor  nodes  cannot  be 
within  the  region  of  the  target  body,  but  the  target  nodes  can  be  inside  or 
outside  the  contactor  body. 

In  the  iterative  solution,  the  conditions  of  sticking  contact, 
sliding  contact,  and  tension  release  at  the  contactor  nodes  for  the 
beginning  of  the  next  iteration  are  determined  from  the  contact 
conditions  of  the  surface  segments  (see  Section  2.6).  The  necessary 
geometric  constraints  to  enforce  the  condition  of  contact  for  the 
incremental  displacements  at  a  generic  contactor  node  are  discussed  in 
this  section.  The  calculation  assumes  that  the  solution  response  is 
known  at  time  t  and  that  (i-1)  iterations  have  been  performed  to 
calculate  the  solution  at  time  t+At. 

Figure  3  shows  how  a  contactor  node  k  has  come  into  contact  with  the 
target  segment  j  formed  by  nodes  A,  B,  C,  and  0.  Defining  point  P  to  be  the 
physical  point  of  contact  of  node  k  in  the  triangle  ABO,  such  that 


tm2p(t-i)  ,  i-i)  .  wt^il  (4) 

where, 

°  current  Qlokal  coordinates  of  node  k  after 
iteration  (i-1)  for  the  equilibrium  configuration 
corresponding  to  time  t+At 

*  rcateHal  overlap  at  contactor  node  k.  The  calculation 
of  overlap  is  such  that  the  vectors  and 

n.T  are  parallel  to  each  other. 
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♦+Atx(i-I), 

t4>At“B(i-|), 

*0 


t+Atx  (i-l), 

“P 


OVERLAP, 


t+At_  (1-1) 


(PARALLEL  TO  nj) 


:smm 


•t+At,  (i-l) 

“A 


W«llhi<whlW  JUki’ktiti‘  \t  t  litif  Jl  *  *^*'*~*  %^4k4M 


wmlsm 


if  / 

f  1  nj  Wigi-  :■& 


TARGET  SEGMENT,  J 


TARGET 

BODY 


Figure  3  Contactor  node  k  in  contact  with  triangle  ABO  of  target 
segment  j  (node  k  has  penetrated  into  the  target  body 
through  triangle  ABO  of  target  segment  j).  Note  that.  In 
general,  the  contactor  node  can  come  into  contact  with  a*;* 
of  the  four  triangles  of  target  segment  j. 


2.2.1  Constraints  on  Surface  Displacements  Due  to  the  Condition  of 
Sticking  Contact  at  Node  k  ~~~~ 

For  a  contactor  node  k  in  sticking  contact,  the  incremental 
displacements  in  iteration  (i)  are  such  that  (see  Fig.  4) 


•  the  material  overlap,  is  eliminated 

•  the  physical  point  of  contact  (point  P  in  Fig.  3)  with  the 
target  segment  remains  unchanged 


Then, 


t+At  (i)  t+At  (i) 

h  "  2p 


(5) 


Subtracting  from  both  sides  of  Eq.  (5), 

t+At  (i)  _  t+At  (1-1)  B  t+At  (i )  _  t+At  (1-1) 


2k 


-P 


X-P 


(6) 


Substituting  Eq.  (4)  Into  Eq.  (6)  and  rearranging,  we  obtain, 

V1’  •  ♦  **y M> 


where, 


(7) 


AUp^  ^  *  incremental  displacement  of  point  P  in  iteration  (1) 

AUj^^  »  incremental  displacement  of  node  k  in  iteration  (i)  . 

Using  Eq.  (3)  and  assuming  an  isoparametric  interpolation,  we  hence  obtain 
the  following  constraint  equation, 

i[4(AuJ<)+^VM)>-<4-(M)  +Y(M))AUA(1) 

“  Y(1'1,AUC(i)  "  1  “  0  * 

(8) 
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TARGET  SEGMENT  j 


CONTACTOR 
NODE  k 


TARGET  SEGMENT  ] 
AFTER  ITERATION 
(i-l) 


MATERIAL  OVERLAP, 

AT  PENETRATED  CONTACTOR  NODE,  k 
IS  PARALLEL  TO  n[ ) 


Figure  4  Condition  of  sticking  contact  at  contactor  node  k 
(  also  see  Fig.  ic  ) 
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where, 


a 


(1-1) 


» 


triangular  area  coordinates  of 
point  P- after  iteration  (i-1) 


incremental  displacements  at  target 
nodes  A,  B,  C,  and  D  respectively, 
in  iteration  (i) 


Equation  (8)  is  the  constraint  of  compatible  surface  displacements  for 
sticking  contact  between  the  contactor  node  k  and  the  target  segment  j. 
Equation  {8}  is  used  in  the  calculation  of  the  incremental  equations  of 
equilibrium  [1], 


2.2.2  Constraint  on  Surface  Displacements  Due  to  the  Condition 
of  Sliding  Contact  at  Node  k 


For  a  contactor  node  k  in  sliding  contact,  the  incremental 
displacements  of  iteration  (i)  are  such  that 

•  the  material  overlap,  t+At^(i-l),  ei iminated 

•  the  physical  point  of  contact  with  the  target  segment  can  change. 


The  area  coordinates  of  the  physical  point  of  contact  after  iteration  (i) 
are  unknown.  So,  assuming  the  amount  of  sliding  in  iteration  (i)  to  be  small 
and  linearizing  about  the  geometry  after  iteration  (1-1),  an  approximate 
constraint  of  compatible  surface  displacements  for  sliding  contact  is  obtained 
as 


(5jT)>y<> 

Substituting  Eq. 


Using  Eq. (3)  and 


-  i  <2jT)>V'>  -  t+4t*p(i-D]  (9) 

(4)  into  Eq.  (9)  and  rearranging,  we  have 

<5jT)TU“p(1)3  =  <n/)V°  *  t+4t*k(M)]  (10) 

assuming  an  isoparametric  interpolation,  we  hence  obtain 
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i  (njT!T[4<iuk(1)  ♦  t+«ik(<-D)  -  (4a'*-1!  ♦  Y<1-I>)auA(,) 

-  (48(1-1)  +Tt1'1))4uB(1)  -T(1'I)4uc(1)-Tt1‘l)MD(1)]  =0-  (11) 

Equation  (11)  represents  a  constraint  equation  on  surface  displacements 
for  sliding  contact  between  the  contactor  node  k  and  the  target  segment  j. 
Equation  (11)  is  used  in  the  calculation  of  the  incremental  equations  of 
equilibrium  [1]. 


2.2.3  Condition  of  Tension  Release  at  Node  k 

For  a  contactor  node  k  which  experiences  tension  release  after 
iteration  (i-1),  the  incremental  nodal  displacements  of  iteration  (i)  are 
independent  of  the  target  segment  displacements. 


Incremental  Equations  of  Equilibrium  for  Contact 


The  governing  equations  of  motion  prior  to  contact  are  derived  using 
the  procedures  described  in  Reference  3.  The  solution  response  of  the 
contactor  and  target  bodies  is  independent  of  each  other  and  the  unknowns  to 
be  calculated  in  each  iteration  are  the  incremental  nodal  point  displacements. 


After  the  occurrence  of  contact,  the  effects  of  constraints  on  incremental 
surface  displacements  (see  Eqs.  (8)  and  (11)),  are  included  in  the  equilibrium 
equations  using  a  Lagrange  multiplier  technique  [1],  The  unknowns  to  be 
calculated  in  each  iteration  are  the  incremental  nodal  point  displacements 
and  the  incremental  contact  forces  (however,  the  calculated  numerical  values 
of  these  incremental  contact  forces  are  not  used  subsequently  in  any  contact 
calculations,  see  Section  4).  The  equilibrium  equations  for  iteration  (i) 
at  time  t+at  are: 
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t+Atjjt'i-l)  B  ysufl^  tangent  stiffness  matrix  Including  material 
and  geometric  nonlinearities  after  iteration  (1-1) 

t+atK  (i-1)  B  Usuaj  tangent  stiffness  matrix  for  the  contactor  body 
after  Iteration  (i-1) 

t+at^  (1-1)  B  Usuaj  ta ngent  stiffness  matrix  for  the  target  body 
after  Iteration  (1-1) 

t+at»  (1-1)  m  Qontact  matrix  to  Include  the  constraints  of 

compatible  surface  displacements  after  iteration 
(i-1)  (see  Section  2.2) 

«  Vector  of  total  applied  external  forces  at  time  t+ot 
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t+Atp(i-l) 


t+Atp  (i-1) 
-c 

t+AtA  (i-1) 


AU 


(i) 


AA 


(i) 


Vector  of  nodal  point  forces  equivalent  to  element 
stresses  after  iteration  (i-1) 

Vector  of  updated  contact  forces  after  iteration 
(i-1)  (see  Section  2.5) 

Vector  of  material  overlaps  at  contactor  nodes  after 
iteration  (i-1)  (see  Section  2.2) 

Vector  of  incremental  nodal  point  displacements  in 
iteration  (i) 

Vector  of  increments  in  contact  forces  in  iteration  (i) 
(see  Section  4). 


In  Eq.  (12),  the  second  row  corresponds  to  the  constraints  on  incremental 
surface  displacements  due  to  contact.  The  increments  in  contact  forces, 

AA^,  enforce  the  constraints  of  contact  for  iteration  (i)  (i.e.,  AA^ 

are  Lagrange  multipliers). 


All  surface  nodes  which  belong  to  the  contact  region  contribute  to  the 

contact  matrices,  »  and  The  sum  total 

effect  of  all  surface  nodes  on  the  contact  matrices  is  obtained  by  summing 
the  individual  contributions  using  the  direct  stiffness  method. 

The  calculation  of  t+AV’*^,  t+AtR,  and  is  performed 

using  the  usual  procedures  (see  Reference  3). 


2.4  Contact  Matrices  for  Sticking  Contact  and  Sliding  Contact 

In  this  section,  the  matrices  (i-1)  t+At^  (i-1)  and  t+At.  (i-1) 

are  given  for  a  generic  region  of  contact  consisting  of  the  contactor  node  k 
and  its  target  segment  of  contact  3  (see  Fig.  3). 
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2.4.1  Condition  of  Sticking  Contact 

The  matrices  for  sticking  contact  are: 

.  Y(i-D, 

*y(M) 

*y<M) 


t+At„  (i-1) 

h 

(15  x  3) 


t^t?  (i-1)  B  _  i(48(i-l)  +  yM-Dj  t+4tj  (i-1) 
**c  -k 

(15  x  1) 

_  iv(i-U  t+At.  (i-1) 
*Y  2j( 


iji-l)  t+At.  (i-1) 


where. 


identity  matrix;  (3  x  3) 


t  &e  (i-l)  a  Contact  force  at  contactor  n<  te  k  after  iteration 

(i-1)  corresponding  to  the  x,y.  and  z  coordinate  axes 
(see  Section  2.5). 


Also,  the  solution  vector  for  iteration  (i)  is 


where  is  the  vector  of  increments  in  the  contact  force  at  contactor 

node  k  in  iteration  (i). 


Note  that  by  substituting  Eqs.  (14),  (15), and  (18)  into  Eq.  (12)  the 
constraint  of  compatible  surface  displacements  for  sticking  contact  is 
generated  (see  Eq,  (8)).  Thus,  to  enforce  sticking  contact,  three  individual 
equations  are  necessary  to  constrain  the  x,y, and  z  incremental  displacements 
of  node  k  to  the  x,y,and  z  incremental  displacements  of  point  P  respectively 
(see  Fig.  4). 


2.4.2  Condition  of  Sliding  Contact 


The  matrices  for  sliding  contact  are: 


t+Aty  (i-1)  = 

h 


(15  x  1) 


“-j 


i  *  S'-")  n/ 


i  n,-T 


1  t(m)  n/ 


1  .I*'1*  SjT 


(19) 


.  J  („  T)T  t+4t.  0-l)j 


(1*  1) 


(20) 


The  vector  t+itR.  (i"U  for  sliding  contact  is  as  given  in  Eq.  (15).  The 


solution  vector  for  iteration  (i)  is 
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AjjW 

(16  X  1) 


V0 


(0 


«  ('> 
k 


where  AA, v  '  is  the  increment  in  the  normal  contact  force  at  contactor  node  k 
in  iteration  (i).  Or, 


AA  W 

-k 


=  Mk<,)JllT 


Note  that  by  substituting  Eqs.  (19),  (20)  and  (21)  into  Eq.  (12),  the 
constraint  of  compatible  surface  displacements  for  sliding  contact  is  generated 
(see  Eq.  (11)).  Thus,  to  enforce  sliding  contact,  one  equation  is  necessary  to 

constrain  the  incremental  displacement  of  node  k  along  the  direction  n.T  to  the 

T  « 

incremental  displacement  of  point  P  along  the  same  direction  n^  (see 
Section  2.2). 


2.5  Evaluation  of  the  Contact  Forces  After  Iteration  (1-1) 


After  each  iteration,  the  generated  contact  forces  at  the  contactor  nodes 
are  updated  such  that  Coulomb's  law  of  friction  is  enforced  in  a  global  sense 
over  each  contactor  segment  [12].  The  contact  forces  at  the  target  nodes  are 
updated  such  that  they  are  statically  equivalent  to  the  updated  forces  at  the 
contactor  nodes.  The  updated  nodal  contact  forces  are  the  elements  of  the 

vector  t+atR  (’“1^  (see  (15)). 
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For  a  generic  contactor  node  k  which  belongs  to  the  region  of  contact, 

the  contact  nodal  point  force  vector  prior  to  updating,  AR.  is 

obtained  as  ~ 


aR  (1-1)  =  t+Atp  (i-1)  _  t+Atf| 
-k  ■— k  -k 


(23) 


where  anc[  are  the  elements  of  vectors  t+Atp(i"l)  and 

t+AtR  corresponding  to  the  degrees  of  freedom  of  node  k  respectively. 

Considering  all  contactor  nodes  belonging  to  the  contact  region,  the 
vector  of  nodal  point  contact  forces  prior  to  updating  is  denoted  as 

AR^1"^  and  the  elements  of  AR^“^  corresponding  to  node  k  are  denoted  as 

45k<M)- 

Physically,  aR^1"^  is  (minus)  the  out-of-balance  force  vector  usually 
encountered  in  nonlinear  analyses  without  contact  conditions.  In  the  Dresence 
of  contact,  when  convergence  is  reached  after  iteration  (i-1),  AR^"*' 
is  equal  to  the  contact  nodal  point  force  vector  at  node  k  (see  Fig.  Id). 

The  updated  contact  forces  after  iteration  (i-1),  t+AtRc^”1^,  are 

calculated  from  the  contactor  surface  nodal  forces,  AR^1’1^,  in  the  following 
three  steps: 


•  Distributed  segment  tractions  are  recovered  on  the  contactor  surface 
such  that  they  are  equivalent  (in  the  virtual  work  sense)  to  the  nodal 

contact  forces,  AR^“^. 


•  The  distributed  contactor  segment  tractions  are  updated  to  satisfy 
Coulomb's  law  of  friction.  The  updated  contactor  surface  nodal  forces 
are  obtained  as  the  consistent  nodal  loads  corresponding  to  the 
updated  segment  tractions. 


•  The  target  surface  updated  nodal  contact  forces  are  obtained  from  the 
contactor  surface  updated  nodal  forces  by  considering  static  equilibrium 
of  the  contact  region  after  iteration  (i-1)  (see  Fig.  3  and  Eq.  (15)). 


2.5.1  Recovery  of  Segment  Tractions  on  the  Contactor  Surface 


i 


IS 

I 


1 


ESS 


The  segment  traction  recovery  calculation  uses  the  following  two 
assumptions: 


•  The  interpolation  of  tractions  over  each  contactor  segment  is 
bilinear. 

•  The  tractions  are  continuous  across  the  segment  boundaries. 


Figure  5  shows  the  distribution  of  segment  tractions  over  a  generic 
contactor  segment  j.  The  consistent  nodal  loads  corresponding  to  the 
distributed  segment  tractions  are  given  by 


4Bj  ■  Sij 


where. 


h 


s 


(24) 


(25) 


k 

nodal  point  values  of  the  segment  tractions  (e.g.  t*  is  th< 
x-component  of  the  segment  traction  at  node  k) 
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m 

t 


Figure  5  Contactor  segment  traction  distribution  (  t*,  t*,  tm,  and  tn 
are  the  values  of  segment  tractions  at  segment“hodal  points 
k,  lt  m,  and  n  respectively  ) 


13 

8 

wit 


arJ 

AR* 

AR, 

X 

y 

AR* 

ar* 

AR. 

X 

y 

ar!" 

arI! 

AR1 

X 

y 

ar; 

AR" 

AR1 

X 

y 

(26) 


=  consistent  nodal  point  forces  corresponding  to  the  segment 

tractions  (e.g.,  AR*  is  the  x-component  of  the  consistent 

nodal  load  at  node  k  due  to  the  distributed  segment  tractions 

over  segment  j  only.  The  total  force  at  node  k  is 

the  sum  of  contributions  from  the  tractions  acting  over  all 
segments  adjoining  node  k). 

G  =  coefficient  matrix  relating  nodal  values  of  segment  tractions 
to  the  corresponding  consistent  nodal  loads. 

The  matrix  G  is  evaluated  by  (2  x  2)  numerical  integration  [3]  as 


where , 


G 

fl 

,:ri, 

JH 

«*» 

hll 

h^ 

2!l3 

Jh 

H 

8 

h2  2 

^23 

^24 

h33 

h34 

symmetric 

h44 

8 


(27) 


(28) 


matrix  of  values  of  the  bilinear  Interpolation 
functions  at  the  (2  x  2)  Gauss  Integration  points 


(29) 


=  diagonal  matrix  of  values  of  the  Jacobian  determinant 
at  the  (2  x  2)  Gauss  integration  points 


hn  =  h22  *  h33  =  h44  =  0.62200847 

hj2  =  ^23  =  ^34  =  ^14  s  0.16666666 
hj3  =  h24  =  0.04465820 


Using  Eq.  (24)  and  summing  the  contributions  from  all  contactor  segments 
belonging  to  the  contact  region,  a  coefficient  matrix  relating  the  nodal  values  of 
segment  tractions  to  the  nodal  contact  forces  (i.e.,  the  contact 
forces  given  by  Eq.  (23))  is  constructed.  A  Gauss  elimination  solution  is 
subsequently  performed  to  obtain  the  unknown  nodal  values  of  the  segment 
tractions. 


2.5.2  Friction  Update  of  Segment  Tractions 


Using  the  recovered  segment  tractions,  the  total  segment  contact  force, 
T^,  is  obtained  from  Eq.  (24)  as: 


Ij  -  C  £T  ij  JT 


(30) 


where 


hU°l 

+ 

h12J2 

+ 

h13J3 

+ 

h14J4 

TijjOj 

+ 

"h22J2 

+ 

"h23°3 

+ 

h24J4 

h13Jl 

+ 

^23J2 

+ 

^33°3 

+ 

"h34J4 

Vi 

+ 

^24J2 

+ 

*34J3 

+ 

W4 
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(31) 


(32) 


Also, 

V  ■  mV  s'l  s' 
itJ  ■  IJ-V 

where  TnJ  and  T^J  are  the  total  normal  and  tangential  segment  contact  forces 
respectively. 

Using  TnJ  and  the  procedure  for  updating  the  segment  tractions  to 
enforce  Coulomb's  law  of  friction  is  as  follows: 

•  Condition  of  Tension  Release 

A  contactor  segment  experiences  tension  release  after  iteration  (i-1)  if 

the  total  normal  segment  contact  force  is  tensile. 

r 

Since  the  segment  normal  vector  n.  points  into  the  continuum  of  the 

J 

contactor  body  (see  fig.  2),  a  tensile  normal  segment  force  acts  in  the 

r 

opposite  direction  of  n j  .  Thus,  segment  tension  release  is  detected  if 
(Ij)T  n.C  <  0  .  (34) 

The  segment  tractions  are  updated  to  zero.  Hence 

lj  =  £  (35) 

Aj  =  0  (36) 

♦ 

where 

A 

t.  =  matrix  of  updated  nodal  point  values  of  segment  tractions. 

A ^  =  matrix  of  consistent  segment  nodal  point  forces 

corresponding  to  the  updated  tractions  over  segment  j. 
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•  Condition  of  Sliding  Contact 


A  contactor  segment  experiences  sliding  contact  after  iteration  (i-1) 
if  the  total  segment  tangential  force  exceeds  the  total  segment 
frictional  capacity,  or 

Ttj  >  Tfj  (37) 

where 

V  ■  I  It  I  (38) 

V  *  »  [(Ij)T n^]  (39) 

=  total  segment  frictional  capacity 

u  *  coefficient  of  friction 

The  tangential  component  of  the  segment  tractions  is  updated  to  be  a 

A 

constant  value  tt  for  the  entire  segment  j  such  that  [1], 


where 


And  also 


J1  +  J2  +  J3  *  J4 
area  of  segment  j 


A  4  A  * 

t,  *  t  J  +  V 

-j  -n  -t 


ij 


(40) 


(41) 


(42) 

(43) 
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where 


tnJ  •  Ctj  n/]  fn//  (44) 

=  nodal  point  values  of  the  normal  component  of  the 
recovered  segment  tractions 

A  i  _  A  A  A  A  T  , 

It  =  C  it  -t  -t  -t  -  (45) 

=  nodal  point  values  of  the  updated  tangential 
segment  traction  (see  Eq.  (40)). 


Note  that  using  Eq.  (40),  the  magnitude  of  the  total  tangential  segment 
force  is  scaled  down  to  equal  the  segment  frictional  capacity.  The 

direction  of  the  updated  tangential  force  is  as  for  the  force  TtJ. 

However,  the  direction  of  the  actual  relative  sliding  of  the  contactor 
segment  j  (with  respect  to  the  target  surface)  in  iteration  (i)  is  not 
enforced  to  be  opposite  to  the  direction  of  the  updated  tangential  force 
(see  Eqs.  (12)  and  (19))* 

For  the  condition  of  contact  with  no  friction,  all  contactor  segments 
which  do  not  experience  tension  release  are  in  sliding  contact. 


•  Condition  of  Sticking  Contact 


A  contactor  segment  experiences  sticking  contact  after  iteration  (1-1) 
if  the  total  segment  tangential  force  is  less  than  (or  equal  to)  the  total 
segment  frictional  capacity,  or 

Ttj<Tfj.  (46) 

The  segment  tractions  satisfy  Coulomb's  law  of  friction  and  thus 


h  ■  (47) 

Aj  »  ARj  ,  (48) 
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For  the  condition  of  contact  with  infinite  friction,  all  contactor 
segments  which  do  not  experience  tension  release  are  in  sticking 
contact. 

By  summing  the  updated  segment  nodal  forces,  the  total  updated  contact 
forces  at  the  contactor  nodes  are  obtained.  For  a  contactor  node  k,  the 

updated  force  is  denoted  as  t+At^(i-l)  (see  £q#  (15)), 

In  general,  the  contact  region  comprises  a  number  of  contactor 
segments  and  solitary  contactor  nodes  (see  Fig.  6).  A  contactor  node  k  is 
denoted  as  a  solitary  node  in  contact  if 


•  The  node  k  is  in  contact 

•  None  of  the  adjoining  contactor  segments  to  node  k  are  in  contact. 

This  condition  occurs  if  for  each  of  the  adjoining  segments,  the 
number  of  segment  nodes  in  contact  is  less  than  four  (see  Section  2.1). 


The  update  of  the  contact  force  at  a  generic  solitary  node  in  contact,  k, 
is  performed  as  follows: 


•  Obtain  the  contact  force  ARjJ1"^  at  the  solitary  node  (see  Eq.  (23)). 

•  At  the  solitary  node,  evaluate  the  normal  vector  to  the  contactor  surface, 

c 

nk  ,  as  the  average  of  the  normal  vectors  of  the  adjoining  surface 
segments . 

L 

•  Calculate  the  normal  contact  force,  T_  ,  and  the  tangential 

k 

contact  force,  ,  as 


[<4£k(M)>  2kC]4C 


0-1)  ,  k 


In 


(49) 

(50) 


•  Evaluate  the  updated  contact  force  at  node  k  in  analogy  to  the 
evaluation  for  a  contact  segment: 


For  tension  release 


t+At.  (i-1) 


»  0 


(51) 
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CONTACTOR  SEGMENTS 
IN  CONTACT  (ALL  NODES 
BELONGING  TO  THE  SEGMENTS 
ARE  IN  CONTACT) 


CONTACT  ( NONE  OF 
THE  ADJOINING  SEGMENTS 
ARE  IN  CONTACT) 


Figure  6  Contact  region  comprised  of  six  contactor  segments 
In  contact  and  two  solitary  nodes  In  contact 
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For  sliding  contact 


t+At.  (1-1) 
For  sticking  contact 


In 


(52) 


t+Atik(1'I)  ■  4Rk(M)‘ 


(53) 


2,6  Evaluation  of  the  Conditions  of  Sticking,  Sliding  and  Tension  Release  at 
the  Contactor  Nodes 


Once  the  conditions  of  the  contactor  segments  have  been  decided  as 
discussed  in  Section  2.5.2,  the  algorithm  determines  the  conditions  of  the 
contactor  surface  nodes  as  shown  in  Table  1. 

For  the  solitary  nodes  in  contact,  the  conditions  at  the  nodes  are  given 
directly  by  the  calculation  of  friction  update. 

A  special  case  arises  when  a  contactor  node  comes  into  contact  with 
the  target  surface  after  iteration  (i-1),  while  the  node  was  not  in  contact  after 
iteration  (1-2).  Its  condition  for  iteration  (i)  is  obtained  ass 


•  sliding  contact  if  the  contact  surfaces  are  frictionless 


•  sticking  contact  if  the  contact  surfaces  are  frictional. 

Using  the  conditions  at  the  contactor  nodes  for  iteration  (1),  the 

matrices  and  ttatAc^"^  are  evaluated  (see  Section  2.4).  For 

iteration  (i),  the  total  number  of  contact  equations  is  equal  to  the  number  of 
sliding  contactor  nodes  plus  three  times  the  number  of  sticking  contactor  nodes. 


III.  SOLUTION  METHOD  FOR  CONTACT-DYNAMICS 

The  solution  of  dynamic  contact  problems  is  obtained  using  the  procedure 
discussed  above  for  static  analysis,  but  now  including  the  effects  of  inertia 
and  damping  forces.  Thus,  the  incremental  equations  of  equilibrium  for 
iteration  (i)  at  time  t+At  are: 
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Table  1  STATE  OF  CONTACTOR  NuDE  AS  DECIDED  EY 
STATES  OF  ADJOINING  CONTACTOR  SEGMENTS 


STATE  OF  ADJOINING  CONTACTOR  SEGMENTS 


STATE  OF 
CONTACTOR  NODE 


One  adjoining  segment 

Other  adjoining  segments 

Sticking 

Sticking 

Sliding 

Tension  release 

Sliding 

Sliding 

Tension  release 

Sticking 


Sliding 


Tension  release 


Tension  release 


Tension 

release 


M  0 
0  0 


AU(1> 

0 

f 

C  0 

0  0 

au^ 

0 

t 

+ 

mJ 

- 

L. 

0  t+Atj^(i-l) 


jt+Atly  (i“l) 
-X 


AU 


(i) 


AX 


(i) 


t+At, 


M  0 
0  0 


(t+At  JJ  ( i — 1 ) 
0 


C  0 
0  0 


(t+At^  ( i  - 1 ) 
0 


t+Atp ( 1 - 1 ) 
0 


t+AtR  (i-1) 


t+At.  (i-1) 

4 


(54) 


where, 


AU 


(1) 


4U(,) 

t+At^i-l) 

t+Atg(i-l) 

M 

C 


=  Vector  of  incremental  velocities  in  iteration  (1) 

=  Vector  of  incremental  accelerations  in  iteration  (i) 

*  Vector  of  velocities  after  Iteration  (i-1) 

a  Vector  of  accelerations  after  iteration  (1-1) 

*  Mass  matrix  of  the  contactor  and  target  bodies 


■  Viscous  damping  matrix  for  the  continua  of  the 
contactor  and  target  bodies  • 


The  evaluation  of  the  contact  matrices  and  is 

performed  as  in  static  analysis. 


For  the  contactor  surface  nodes  belonging  to  the  region  of  contact,  the 
contact  forces  after  iteration  (1)  (prior  to  updating)  are  obtained  as 


AR 


(i-1)  a  t+Atp ( 1 — 1 )  +  M  t+Atj|(  1-1)  +  c  t+Atfl(i-l)  _  t+AtR 


(55) 


This  vector  is  used  as  described  in  Section  2.5  to  evaluate  the  updated 


nodal  point  contact  forces 
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IV.  ITERATION  STRATEGY  AND  CONVERGENCE  CRITERIA 

For  iteration  (i)  at  time  t+At,  the  sequence  of  calculations  performed 
in  the  contact  solution  is  as  follows: 

•  Use  the  nodal  forces  aR^1"^  to  recover  segment  tractions  on  the 
contactor  surface.  Update  the  segment  tractions  and  the  forces  acting 
at  the  solitary  nodes  in  contact  to  satisfy  Coulomb's  law  of 

friction  and  calculate  the  updated  contact  surface  forces  t+AtR 
(see  Section  2.5).  ~c 

•  Use  the  current  geometry  after  iteration  (i-1)  to  determine  the 
target  segment  of  contact  for  each  contactor  node  belonging  to  the 
contact  region.  Also  evaluate  the  material  overlap  andthe  area 
coordinates  of  the  physical  point  of  contact. 

•  Detect  if  any  new  contactor  nodes  have  come  into  contact  after 
iteration  (i-1). 

•  Evaluate  the  states  of  the  contactor  nodes  for  iteration  (i)  (see 
Section  2.6). 

•  Evaluate  the  contact  matrices  /i_1)  and  and 

assemble  all  matrices  as  given  in  Eq.  (12).  Then  solve  Eq.  (12) 
to  obtain  the  unknown  quantities  AU^  and  aa^. 

Note  that  the  increments  in  contact  forces,  aa^\  obtained  from  the 
solution  of  Eq.  (12)  are  not  used  subsequently  in  any  calculation.  The 
total  contact  forces  (prior  to  updating)  after  each  iteration  are  simply 
obtained  as  given  by  Eq.  (23). 

Convergence  of  solution  is  accepted  after  iteration  (i)  if  the  following 
criteria  are  simultaneously  satisfied: 


•  Energy  convergence  criterion 


AU{i)T[t+AtR  -  t+Atp(M)  .  M  t+At(j(i-i)  .  Ct+Aty(i-1)  +  t+AtR  (i-l). 

. .  .  «**—  J 

~ - - - - - - “ - <ET0L 

AU(1)  j-t+AtR  „  tp  .  H  t+Aty(o)  .  c  UAt^(o)  +  ^  •, 

(56) 
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where  ETOL  is  an  energy  convergence  tolerance. 


•  Force  convergence  criterion 


JJt+Atpj  _  t+Atp(i-l)  _  M  t+Atfl(i-l)  _  c  t+At^i-l)  +  t+AtR  (i-1) 

RNORM 


S  RTOL 
(57) 


where  RNORM  is  a  reference  force  for  convergence  [7]  and  RTOL  is  a  force 
convergence  tolerance.  In  the  evaluation  of  the  Euclidean  norm  in  Eq.  (57), 
only  the  translational  degrees  of  freedom  are  considered. 


For  the  rotational  degrees  of  freedom,  a  moment  convergence  criterion 
is  specified  in  analogy  to  the  force  convergence  criterion. 


•  Contact  Force  Convergence  Criterion 


4R'1-1)  .  arCI-2)  I, 


AR 


(1-1) 


<  RCTOL 


+  RCONSM 


(58) 


where  RCTOL  is  a  contact  force  convergence  tolerance  and  RCONSM  is  a 
small  number.  Note  that  RCONSM  makes  the  denominator  of  Eq.  (58) 

nonzero  when  j|  AR^’^  jjg  ■  0  (i.e.,  no  contact  conditions  exist 

during  iteration  (i-1)). 

V.  TIME  INTEGRATION  OF  THE  EQUILIBRIUM  EQUATIONS  FOR  DYNAMIC  CONTACT 

A  valid  dynamic  contact  solution  must  satisfy  the  following  two  criteria: 

•  The  total  energy  of  the  system  is  conserved 

•  The  impulse-momentum  relationship  is  satisfied  for  each  body. 


Using  the  Newmark  method  of  implicit  time  integration,  the 
velocity  vectors  after  iteration  (i-i)  at  time  t+At  are  [3] 


acceleration  and 


“s""1  ■  slip  ‘“V’-'-'.J  -  -dr'i  •  < i* 


.  'u  *  [(l.i)'K  .  (  PPyP-ll]  „  ,60, 

where, 

a, 6  -  Newmark  parameters 
t  t*  t” 

U,  U,  U  =  displacement,  velocity,  and  acceleration  vectors  at 
time  t  respectively 

t+At,,(i-l)  i. 

U  =  displacement  vector  after  iteration  (i-1)  at  time  t+At  , 

Also,  the  incremental  acceleration  and  velocity  vectors  in  iteration  (i)  are 
given  by  ' 


4u<(> 


AU«> 

(61) 

*U«>. 

(62) 

Equations  (59)  to  (62)  are  used  in  the  incremental  equations  of  equilibrium 
for  dynamic  contact  (see  Eq.  (54)).  Since  the  calculated  response  depends 
on  the  values  of  the  Newmark  parameters  a  and  6,  the  objective  must  be  to 
choose  optimum  parameters  for  the  solution. 

In  search  for  suitable  a  and  4  parameters,  the  solutions  to  some  simple 
numerical  experiments  of  two  unequal  colliding  point  masses  were  considered. 


The  experiments  were  repeated  for  various  combinations  of  Impact  velocities. 
The  study  showed  that  the  criteria  of  energy  and  momentum  balance 
are  satisfied  when  and  6=i  are  used.  With  other  values  of  o  and  6, 
the  energy  and  momentum  equations  are  not  necessarily  satisfied. 

Another  numerical  experiment  Involves  the  longitudinal  impact  of  two 
identical  bars  which  are  moving  towards  each  other  with  a  constant  velocity 
(see  Fig.  7a, b).  The  bar  material  is  linear  elastic  and  the  contact 
surfaces  are  frictionless.  For  the  solution  obtained  using  a=i  and  6=i, 
the  total  energy  of  the  bar  (strain  energy  +  kinetic  energy)  is  conserved 
throughout  the  solution  and  the  impulse-momentum  relationship  is  satisfied. 
Good  agreement  is  obtained  between  the  numerical  solution  and  the  analytical 
solution  [4]  for  the  stress  generated  due  to  Impact  (see  Fig.  7c). 

In  addition,  the  integration  scheme  with  <x-i  and  has  the 
following  desirable  characteristics  (considering  linear  analysis): 


•  The  method  is  unconditionally  stable  since  the 
condition  o  6  i  (6+i)2  is  satisfied. 


•  The  scheme  gives  no  amplitude  decay,  since  6  =  }. 


•  The  percentage  period  elongation  is  reasonably  small. 


VI.  NUMERICAL  SOLUTIONS 


The  algorithm  presented  in  the  previous  sections  was  implemented  in 
the  computer  program  ADINA  [7]  and  the  results  of  some  sample  analyses  are 
presented  in  this  section.  In  these  analyses,  the  primary  objective  was  to 
study  the  performance  of  the  algorithm  under  various  conditions  of  contact. 


6.1  Analysis  of  Ax i symmetric  Hert2  Contact  Problems 


A  sphere  of  radius  R*5  is  considered  to  come  into  contact  with  a  flat 
rigid  surface  (see  Fig,  8).  The  analysis  Is  performed  for  the  following 
two  conditions: 


•  Static  analysis  of  contact  when  the  sphere  is  subjected  to  an  externally 
applied  total  force  P*-523k  (see  Fig.  8a) 
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Figure  7  Analysis  of  longitudinal  impact  of  identical  bars 


b)  Finite  element  model 


c)  Impact  stress  In  element  number  1.  Each  solution  point  is  the  average  element 

stress  over  a  time  interval  0.0125  seconds.  The  time  step  used  is  At=0. 00025  seconds. 


•  Dynamic  analysis  of  contact  when  the  sphere  impacts  the  rigid 

surface.  Prior  to  impact,  the  sphere  has  an  initial  velocity  V=-3k. 


6.1.1  Finite  Element  Discretization 


Due  to  symmetry  of  the  sphere  geometry  and  the  applied  loading,  only  a 
thin  wedge  from  the  sphere  continuum  is  discretized  (see  Fig.  8b, c).  The  wedge 
is  bounded  by  its  two  semicircular  sides  and  the  enclosed  wedge  angle. is 
1  degree.  The  wedge  is  discretized  using  eight  node  3-D  solid  elements. 

For  each  node  on  Face  1  of  the  wedge,  the  degree-of-freedom  normal 

to  Face  1  ( i . e . ,  along  the  skew  coordinate  direction  y^j  is  deleted.  Similarly, 

for  each  node  on  Face  2  of  the  wedge,  the  degree-of-freedom  normal  to  Face  2 

(i.e.,  along  the  skew  coordinate  direction  y^)  is  deleted. 

The  contactor  surface  is  defined  over  the  wedge  boundary  (see  Fig.  8b).  The 
target  surface  is  defined  to  be  the  flat  rigid  surface,  which  is  modeled  by 
specifying  target  nodes  with  no  degrees  of  freedom.  In  an  X— Y  plane  view,  the 
target  surface  is  triangular  in  shape. 


6.1.2  Static  Analysis 


The  wedge  is  subjected  to  a  uniformly  distributed  body  force  along  the 
negative  z  direction  such  that  the  resultant  of  the  body  force  corresponds 
to  P=-523k  acting  at  the  center  of  the  sphere  (i.e.,  the  total  force  applied 
onto  the  wedge  equals  -( 523/360) k). 

Figure  9  shows  the  calculated  contact  tractions  and  a  comparison  with  the 
Hertz  analytical  solution  [5]. 


6,1.3  Dynamic  Analysis 


For  all  nodes  belonging  to  the  wedge,  an  initial  velocity  V*-3k  was 
assigned.  The  time  integration  of  the  dynamic  response  was  performeH  using 
the  Newmark  parameters  and  6=i  (see  Section  5).  The  time  step  was  equal 
to  At°0.01  sec. 

Figure  10  shows  the  calculated  contact  tractions  and  a  comparison  with  the 
Hertz  analytical  solution  [5]. 


6.2  Analysis  of  Compressed  Spheres  Subjected  to  a  Torsional  Moment 


Figure  11  shows  two  identical  spheres  which  are  first  compressed  onto 
each  other  (P--222k)  and  then  subjected  to  a  twisting  moment(Ma52k). 

The  coefficient  of  friction  is  very  large  (infinite). 
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i  sTi* 


c)  Finite  element  discretization  of  the  wedge  (  Xj— Zj  plane  view  ). 
One  3-D  element  is  used  across  the  wedge  thickness  (  see  Fig.  8b  ) 

Figure  8  Continued 
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Solution  to  the  Hertz  axisymmetric  contact  problem  (statics) 


TOTAL 

LOAD,  P  *-222  k 


Figure  11  Analysis  of  compressed  spheres  subjected  to  a  torsional 
moment 
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Due  to  the  symmetry  in  geometry  and  the  skew  symmetry  of  the  applied 
loading,  we  can  model  the  problem  by  considering  a  single  sphere  first 
compressed  onto  a  flat  rigid  surface  and  then  subjected  to  a  twisting  moment. 
This  sphere  is  modeled  as  a  thin  wedge,  as  in  the  solutions  of  Section  6.1. 

The  finite  element  model  of  the  wedge  is  as  shown  in  Figs.  8b, c.  To 
enforce  the  condition  of  skew-symmetry,  the  x^y^,  arid  Zj  displacements  of 

the  nodes  on  Face  1  are  constrained  to  be  respectively  equal  to  the 
*2»y2’  and  z2  displacements  of  the  corresponding  nodes  on  Face  2. 

The  wedge  is  first  subjected  to  a  uniformly  distributed  body  force  along 
the  negative  z  direction  such  that  the  resultant  of  the  body  force  corresponds 
to  P=-222k  acting  at  the  center  of  the  sphere.  As  a  result,  a  region  of 
contact  between  the  sphere  and  the  rigid  surface  is  established  (radius 
of  contacts. 9). 

To  apply  the  twisting  moment,  the  wedge  is  subjected  to  a  distributed 
lateral  body  force  such  that  at  any  point  (x,y,z)  within  the  wedge,  the  body 
force  per  unit  volume,  o,  is 


b  =  C  r  i  (63) 

-0 

where. 


C  =  constant 


a  distance  of  point  (x,y,z)  from  the  Z  axis 


i 

“Q 


(65) 


«  unit  vector  orthogonal  to  x^  +  yj. 

For  00.01,  the  applied  lateral  body  force  corresponds  to  a  total  twisting 
moment  M“S2k  applied  to  the  sphere. 

Figure  12  shows  the  calculated  shearing  traction  in  the  contact  region 
and  a  comparison  with  the  analytical  solution  [6]. 


/namic  Analysis  of  Frictional  Flidim 
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Figure  13  shows  a  mass,  m=0.2,  attached  to  two  identical  springs  which 
are  anchored  to  a  flat  rigid  surface.  The  rigid  surface  lies  in  the  X— Y 
plane  and  the  mass  rests  on  the  surface.  Thecoefficient  of  friction  between 
the  mass  and  tne  surface  is  0.15.  Both  springs  have  the  same  linear 
force-deflection  relationship  and  geometric  nonlinearities  are  included  in  the 
analysis. 


RIGID  SURFACE 
/tsQ.15 


m*0.2 


L 


SPRINGS 
E  *  1000.0 
A  *1.0 


ACCELERATION 
DUE  TO  GRAVITY  «HOk 


a)  Problem  considered 


Figure  13  Dynamic  analysis  of  frictional  sliding  of  a  point  mass 


The  following  two  cases  are  considered  to  study  effects  of  friction  on  the 
vibrations  of  the  mass-spring  system: 


•  The  point  mass  is  given  an  initial  velocity  of  -l.Oi.  and  is  then  released. 
The  movement  of  the  mass  takes  place  along  the  x-direction  and  the  motion 
is  resisted  by  the  developed  frictional  force  (when  the  mass  is  moving 
along  the  positive  x  direction,  the  frictional  force  acting  on  the 

mass  equals  -0.3(H). 

•  The  point  mass  is  given  an  initial  displacement  of  -2.0j  and  is  then 
released.  The  movement  of  the  mass  takes  place  along  tne  y  direction  and 
the  motion  is  resisted  by  the  developed  frictional  force  (when  the  mass 
is  moving  along  the  positive  y  direction,  the  frictional  force  acting 

on  the  mass  equals  -0.30^) . 

Target  Surface  The  flat  rigid  surface  is  chosen  to  be  the  target  surface  and 
is  modeled  using  four  nodes  each  with  no  degrees  of  freedom. 

Contactor  Surface  The  mass-spring  system  is  chosen  to  be  the  contactor  body. 

An  auxiliary  node  (fixed  in  space)  is  used  to  define  the  contactor  surface 
consisting  of  two  contactor  segments  (see  Figure  13b).  The  auxiliary  node 
lies  outside  the  region  of  the  rigid  target  surface  and  thus  is  not  in 
contact  throughout  the  analysis.  As  a  result,  nodes  1,2,  and  3  are  treated 
as  solitary  nodes  in  contact.  The  direction  of  the  normal  vector  at  the 
solitary  nodes  is  obtained  from  the  geometry  of  the  contactor  segments. 

Figure  14  shows  the  numerical  results  for  the  two  cases  considered. 

The  obtained  solutions  satisfy  the  criteria  of  energy  and  momentum  balance 
and,  for  each  case,  the  work  done  by  the  frictional  force  equals  the 
difference  between  the  initial  energy  and  the  final  energy  of  the  system. 


VII.  CONCLUDING  REMARKS 

An  algorithm  for  the  solution  of  static  and  dynamic  three-dimensional 
contact  problems  has  been  presented.  The  solution  procedure  uses  a  Lagrange 
multiplier  technique  to  incrementally  impose  the  constraints  of  compatible 
surface  displacements  due  to  contact.  The  contact  forces  are  evaluated  from 
distributed  tractions  that  act  on  the  contacting  surfaces.  These  tractions  are 
evaluated  from  the  nodal  point  forces,  which  correspond  to  the  internal  element 
stresses  and  the  externally  applied  loading,  and  the  frictional  conditions 
based  on  Coulomb's  law.  Some  solution  results  obtained  using  the  algorithm 
have  been  presented  to  demonstrate  the  applicability  and  performance  of  the 
solution  procedure. 
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b)  Finite  element  model  used.  All  nodes  except  node  2  have 
all  degrees  of  freedom  deleted. 


Figure  13  Continued 


Using  the  algorithm  already  a  wide  variety  of  contact  problems 
can  be  analyzed.  However,  there  is  an  ever  increasing  need  to  solve  more 
complex  contact  problems  and  additional  research  should  be  performed  to 
broaden  the  applicability  and  effectiveness  of  the  algorithm.  This  research 
should  focus,  for  example,  on  the  following  items: 

•  Line  searching  in  equilibrium  iterations  in  the  presence  of  contact. 

•  Development  of  the  algorithm  for  explicit  time  integration. 

•  Improvement  and  evaluation  of  accuracy  and  effectiveness  of  implicit  and 
explicit  time  integration  schemes  for  dynamic  contact. 

•  Modeling  of  contact  surfaces  using  higher  order  surface  segments. 

•  Evaluation  of  various  friction  laws  and  corresponding  implementation 
for  improved  modeling  of  interface  conditions  and  solution 
effectiveness. 

These  studies  would  be  very  useful  to  gain  a  deeper  understanding  of 
contact  phenomena. 
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